#!pip install scikit-optimize
#!pip install wbgapi
#!pip install xgboost

# Importing Data and Libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import datetime as dt

import wbgapi as wb

import xgboost as xgb

from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import TimeSeriesSplit

np.random.seed(42)
rs = 42

# Importing Data from World Bank

# Specify the indicators to be imported

'''
The Choice of the Following indicators was selected based on the Following Papers:
Qureshi, S., Chu, B. M., & Demers, F. S. (2020). Forecasting canadian gdp growth using xgboost (No. 20-14). Carleton University, Department of Economics.
Adewale, M. D., Ebem, D. U., Awodele, O., Sambo-Magaji, A., Aggrey, E. M., Okechalu, E. A., ... & Oluyide, O. P. (2024). Predicting gross domestic product using the ensemble machine learning method. Systems and Soft Computing, 6, 200132.
'''

indicators = {
    "FP.CPI.TOTL.ZG": "Inflation (CPI annual %)",
    "SL.UEM.TOTL.ZS": "Unemployment rate (%)",
    "NE.CON.PRVT.KD": "Household consumption (constant US$)",
    "NE.EXP.GNFS.KD": "Total exports (constant US$)",
    "NE.IMP.GNFS.KD": "Total imports (constant US$)",
    "BX.KLT.DINV.CD.WD": "FDI net inflows (current US$)",
    "FM.LBL.MQMY.CN": "Money supply M2 (local currency)",
    "CM.MKT.LCAP.CD": "Stock market capitalization (current US$)",
    "SP.POP.TOTL": "Total population",
    "SL.TLF.TOTL.IN": "Total labor force",
    "SP.URB.TOTL": "Urban population",
    "NV.AGR.TOTL.KD": "Agriculture value-added (constant US$)",
    "NV.IND.TOTL.KD": "Industry value-added (constant US$)",
    "NV.SRV.TOTL.KD": "Services value-added (constant US$)",
    "FP.CPI.TOTL": "Consumer Price Index",
    "NE.CON.GOVT.ZS": "Government consumption (real, % of GDP)",
    "NE.CON.PRVT.ZS": "Private consumption (real, % of GDP)",
    "BN.CAB.XOKA.GD.ZS": "Current account balance (% of GDP)",
    "SL.EMP.TOTL.SP.ZS": "Employment (% of total population)",
    "PA.NUS.FCRF": "Official exchange rate (LCU per US$, period average)",
    "NY.GDP.DEFL.KD.ZG": "GDP deflator",
    "NY.GDP.MKTP.KD.ZG": "Real GDP growth (annual %)",
    "NE.IMP.GNFS.ZS": "Imports (real, % of GDP)",
    "NV.IND.TOTL.ZS": "Industrial production index",
    "NE.GDI.FTOT.ZS": "Fixed investment (total, real)",
    "FI.RES.TOTL.CD": "Foreign exchange reserves (current US$)",
    'CM.MKT.LCAP.GD.ZS': 'Market capitalization of listed domestic companies (% of GDP)'
}

# OECD Countries
countries = ["AUS", "AUT", "BEL", "CAN", "CZE", "DNK", "FIN", "FRA", "DEU", "GRC", "HUN", "ISL", "IRL", "ITA", "JPN", "KOR", "LUX", "MEX", "NLD", "NZL", "NOR", "POL", "PRT", "SVK", "ESP", "SWE", "CHE", "TUR", "GBR", "USA"]

# Importing data
raw_data = wb.data.DataFrame(indicators, countries, time=range(1991, 2022)) # From 1991 to 2022
data = raw_data.copy()

# Transforming the data
data = data.stack().unstack(level=1).reset_index().rename(columns={'level_1': 'Year'})

# Converting the Year column to datetime
data['Year'] = data['Year'].str.replace(r'[YR]', '', regex=True)
data['Year'] = pd.to_datetime(data['Year'], format='%Y')
data = data.set_index('Year')
data_no_dummies_raw = data.copy() # Creating a copy of the data without dummies
data = pd.get_dummies(data) # Creating Dummy Variables

# Printing number of rows in the dataset
print(f"Number of Rows in the Dataset: {len(data)}")

# Printing the first 5 rows of the data
print(data.tail())

Number of Rows in the Dataset: 930
            BN.CAB.XOKA.GD.ZS  BX.KLT.DINV.CD.WD  CM.MKT.LCAP.CD  \
Year                                                               
2017-01-01          -1.874424       3.808230e+11    3.212070e+13   
2018-01-01          -2.129357       2.147150e+11    3.043631e+13   
2019-01-01          -2.050828       3.159840e+11    3.408573e+13   
2020-01-01          -2.815417       1.370660e+11    4.156986e+13   
2021-01-01          -3.665292       4.758050e+11    4.854854e+13   

            CM.MKT.LCAP.GD.ZS  FI.RES.TOTL.CD  FP.CPI.TOTL  FP.CPI.TOTL.ZG  \
Year                                                                         
2017-01-01         163.780010    4.512853e+11   112.411557        2.130110   
2018-01-01         147.344853    4.499071e+11   115.157303        2.442583   
2019-01-01         158.243992    5.167006e+11   117.244195        1.812210   
2020-01-01         194.669183    6.283697e+11   118.690502        1.233584   
2021-01-01         205.009026    7.161523e+11   124.266414        4.697859   

            NE.CON.GOVT.ZS  NE.CON.PRVT.KD  NE.CON.PRVT.ZS  ...  economy_MEX  \
Year                                                        ...                
2017-01-01       13.847134    1.293167e+13       67.767473  ...        False   
2018-01-01       13.844213    1.328613e+13       67.457852  ...        False   
2019-01-01       13.983559    1.357142e+13       67.026718  ...        False   
2020-01-01       14.895745    1.322755e+13       66.617903  ...        False   
2021-01-01       14.253033    1.438786e+13       68.045389  ...        False   

            economy_NLD  economy_NOR  economy_NZL  economy_POL  economy_PRT  \
Year                                                                          
2017-01-01        False        False        False        False        False   
2018-01-01        False        False        False        False        False   
2019-01-01        False        False        False        False        False   
2020-01-01        False        False        False        False        False   
2021-01-01        False        False        False        False        False   

            economy_SVK  economy_SWE  economy_TUR  economy_USA  
Year                                                            
2017-01-01        False        False        False         True  
2018-01-01        False        False        False         True  
2019-01-01        False        False        False         True  
2020-01-01        False        False        False         True  
2021-01-01        False        False        False         True  

[5 rows x 56 columns]

data.info()

<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 930 entries, 1991-01-01 to 2021-01-01
Data columns (total 56 columns):
 #   Column             Non-Null Count  Dtype  
---  ------             --------------  -----  
 0   BN.CAB.XOKA.GD.ZS  864 non-null    float64
 1   BX.KLT.DINV.CD.WD  915 non-null    float64
 2   CM.MKT.LCAP.CD     755 non-null    float64
 3   CM.MKT.LCAP.GD.ZS  755 non-null    float64
 4   FI.RES.TOTL.CD     925 non-null    float64
 5   FP.CPI.TOTL        930 non-null    float64
 6   FP.CPI.TOTL.ZG     928 non-null    float64
 7   NE.CON.GOVT.ZS     926 non-null    float64
 8   NE.CON.PRVT.KD     924 non-null    float64
 9   NE.CON.PRVT.ZS     926 non-null    float64
 10  NE.EXP.GNFS.KD     925 non-null    float64
 11  NE.GDI.FTOT.ZS     926 non-null    float64
 12  NE.IMP.GNFS.KD     925 non-null    float64
 13  NE.IMP.GNFS.ZS     926 non-null    float64
 14  NV.AGR.TOTL.KD     849 non-null    float64
 15  NV.IND.TOTL.KD     849 non-null    float64
 16  NV.IND.TOTL.ZS     872 non-null    float64
 17  NV.SRV.TOTL.KD     849 non-null    float64
 18  NY.GDP.DEFL.KD.ZG  930 non-null    float64
 19  NY.GDP.MKTP.KD.ZG  930 non-null    float64
 20  PA.NUS.FCRF        926 non-null    float64
 21  SL.EMP.TOTL.SP.ZS  930 non-null    float64
 22  SL.TLF.TOTL.IN     930 non-null    float64
 23  SL.UEM.TOTL.ZS     930 non-null    float64
 24  SP.POP.TOTL        930 non-null    float64
 25  SP.URB.TOTL        930 non-null    float64
 26  economy_AUS        930 non-null    bool   
 27  economy_AUT        930 non-null    bool   
 28  economy_BEL        930 non-null    bool   
 29  economy_CAN        930 non-null    bool   
 30  economy_CHE        930 non-null    bool   
 31  economy_CZE        930 non-null    bool   
 32  economy_DEU        930 non-null    bool   
 33  economy_DNK        930 non-null    bool   
 34  economy_ESP        930 non-null    bool   
 35  economy_FIN        930 non-null    bool   
 36  economy_FRA        930 non-null    bool   
 37  economy_GBR        930 non-null    bool   
 38  economy_GRC        930 non-null    bool   
 39  economy_HUN        930 non-null    bool   
 40  economy_IRL        930 non-null    bool   
 41  economy_ISL        930 non-null    bool   
 42  economy_ITA        930 non-null    bool   
 43  economy_JPN        930 non-null    bool   
 44  economy_KOR        930 non-null    bool   
 45  economy_LUX        930 non-null    bool   
 46  economy_MEX        930 non-null    bool   
 47  economy_NLD        930 non-null    bool   
 48  economy_NOR        930 non-null    bool   
 49  economy_NZL        930 non-null    bool   
 50  economy_POL        930 non-null    bool   
 51  economy_PRT        930 non-null    bool   
 52  economy_SVK        930 non-null    bool   
 53  economy_SWE        930 non-null    bool   
 54  economy_TUR        930 non-null    bool   
 55  economy_USA        930 non-null    bool   
dtypes: bool(30), float64(26)
memory usage: 223.4 KB

# Define features (X) and target variable (y)
X = data.drop(columns=['NY.GDP.MKTP.KD.ZG'])
y = data['NY.GDP.MKTP.KD.ZG']

# Split Data into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=rs)

# Initialize and train the XGBoost model
model = xgb.XGBRegressor(objective='reg:squarederror', random_state=rs) # Use appropriate objective function
model.fit(X_train, y_train)

# Make predictions on the test set
predictions = model.predict(X_test)

# Evaluate the model (example: using mean squared error)
mse = mean_squared_error(y_test, predictions)
mape = mean_absolute_percentage_error(y_test, predictions)

print("Printing Metrics for Basic XGBoost Model Predicting GDP Growth")
print(f"Mean Squared Error: {mse}")
print(f"Mean Absolute Percentage Error: {mape}")

Printing Metrics for Basic XGBoost Model Predicting GDP Growth
Mean Squared Error: 6.173963391360166
Mean Absolute Percentage Error: 1.4760121339754226

# Pre-Proccessing
data_no_dummies = data_no_dummies_raw.copy()

# Create a lagged variable for GDP Growth
gdp_growth_lag = data_no_dummies.groupby('economy')['NY.GDP.MKTP.KD.ZG'].shift(-1)

# Add the lagged variable to the data
data_no_dummies['gdp_growth_lag'] = gdp_growth_lag
data_no_dummies.dropna(subset=['gdp_growth_lag', 'NY.GDP.MKTP.KD.ZG'],inplace=True) # Dropping Last Observation for each Economy

data_lag = pd.get_dummies(data_no_dummies) # Creating Dummy Variables

# Slip Data into Training and Testing Sets
train = data_lag.loc[data_lag.index <= '2013']
test = data_lag.loc[data_lag.index > '2013']

X_train = train.drop(columns=['gdp_growth_lag'])
y_train = train['gdp_growth_lag']

X_test = test.drop(columns=['gdp_growth_lag'])
y_test = test['gdp_growth_lag']

# Running a hyper parameter sensitivity analysis on Max Depth.

hyper_param_analysis = {}

for i in range(1,30):
  model_max_d = xgb.XGBRegressor(random_state=rs, max_depth=i)
  model_max_d.fit(X_train, y_train)
  predicts = model_max_d.predict(X_test)

  mse = mean_squared_error(y_test, predicts)
  hyper_param_analysis[i] = mse

hyper_param_analysis_df = pd.DataFrame(list(hyper_param_analysis.items()), columns=['max_depth', 'mse'])
hyper_param_analysis_df = hyper_param_analysis_df.sort_values(by='max_depth', ascending=True)

print(hyper_param_analysis_df)

    max_depth        mse
0           1  14.077543
1           2  13.153633
2           3  13.361508
3           4  14.221084
4           5  14.572756
5           6  13.764674
6           7  14.399877
7           8  14.017956
8           9  14.916071
9          10  15.632574
10         11  15.466175
11         12  15.062584
12         13  16.105452
13         14  16.453865
14         15  15.939274
15         16  15.179350
16         17  15.319656
17         18  15.459264
18         19  15.298920
19         20  15.185092
20         21  15.217684
21         22  15.452439
22         23  15.458359
23         24  15.351249
24         25  15.382546
25         26  15.397722
26         27  15.383109
27         28  15.375365
28         29  15.387533

plt.figure(figsize=(12, 8))

# Plotting the Hyperparameter Sensitivity Analysis
sns.lineplot(data=hyper_param_analysis_df, x='max_depth', y='mse')

# Add Points to the Line Plot
plt.scatter(hyper_param_analysis_df['max_depth'], hyper_param_analysis_df['mse'], color='orange')
plt.xlabel('Max Depth')
plt.ylabel('Mean Squared Error')
plt.title('Hyperparameter (Max Depth) Sensitivity Analysis')
plt.grid(True)

# Add Horizontal Line For Best Max Depth
best_max_depth = hyper_param_analysis_df.loc[hyper_param_analysis_df['mse'].idxmin(), 'max_depth']
plt.axhline(y=hyper_param_analysis_df['mse'].min(), color='r', linestyle='--', label=f'Best Max Depth: {best_max_depth} \nMSE: {round(hyper_param_analysis_df["mse"].min(),2)}')

# Add Max Depth Labels and prevent overlapping
for i in range(len(hyper_param_analysis_df)):
    if i < 15: 
        plt.text(hyper_param_analysis_df['max_depth'].iloc[i] + 0.5, hyper_param_analysis_df['mse'].iloc[i] + 0.01, 
                str(hyper_param_analysis_df['max_depth'].iloc[i]), 
                horizontalalignment='left', size='medium', color='black')
    elif i % 2 == 0:
        plt.text(hyper_param_analysis_df['max_depth'].iloc[i] + 0.5, hyper_param_analysis_df['mse'].iloc[i] + 0.01, 
                str(hyper_param_analysis_df['max_depth'].iloc[i]), 
                horizontalalignment='right', size='medium', color='black')


plt.legend(loc='upper right')
plt.show()

# Create a lagged variable for GDP Growth
data_no_dummies = data_no_dummies_raw.copy()

# Create a lagged variable for GDP Growth
gdp_growth_lag = data_no_dummies.groupby('economy')['NY.GDP.MKTP.KD.ZG'].shift(-1)

# Add the lagged variable to the data
data_no_dummies['gdp_growth_lag'] = gdp_growth_lag
data_no_dummies.dropna(subset=['gdp_growth_lag', 'NY.GDP.MKTP.KD.ZG'],inplace=True) # Dropping Last Observation for each Economy

data_lag = pd.get_dummies(data_no_dummies) # Creating Dummy Variables

# Slip Data into Training and Testing Sets
train = data_lag.loc[data_lag.index <= '2013']
test = data_lag.loc[data_lag.index > '2013']

X_train = train.drop(columns=['gdp_growth_lag'])
y_train = train['gdp_growth_lag']

X_test = test.drop(columns=['gdp_growth_lag'])
y_test = test['gdp_growth_lag']


# Create a TimeSeriesSplit object for time series cross-validation
tscv = TimeSeriesSplit(n_splits=5)

# Hyperparameter Tuning
from skopt import BayesSearchCV
from skopt.space import Real, Integer

# Define the hyperparameter grid
param_space = {
    'n_estimators': Integer(100, 5000),
    'max_depth': Integer(3, 15),
    'learning_rate': Real(1e-4, 1e-1, prior='log-uniform'),
    'subsample': Real(0.1, 1),
    'colsample_bytree': Real(0.2, 1),
    'gamma': Real(0, 7),
    'min_child_weight': Integer(1, 10),
    'reg_alpha': Real(1e-8, 10.0, prior='log-uniform'),  # L1 regularization
    'reg_lambda': Real(1e-8, 10.0, prior='log-uniform') # L2 regularization
}

# Initialize the XGBoost model
model_proj = xgb.XGBRegressor(random_state=rs)

# Training the Model with Time Series Cross Validation
bayes_search = BayesSearchCV(estimator=model_proj, search_spaces=param_space, cv=tscv, scoring='neg_mean_squared_error')
bayes_search.fit(X_train, y_train)

# Retrieve the best model and best hyperparameters
best_model_proj = bayes_search.best_estimator_
best_params = bayes_search.best_params_

# Make predictions on the test set
predictions = best_model_proj.predict(X_test)

# Evaluate the model (example: using mean squared error)
mse = mean_squared_error(y_test, predictions)
mape = mean_absolute_percentage_error(y_test, predictions)

print("Printing Metrics for XGBoost Model Forcasting GDP Growth")
print(f"Mean Squared Error: {mse}")
print(f"Mean Absolute Percentage Error: {mape}")

Printing Metrics for XGBoost Model Forcasting GDP Growth
Mean Squared Error: 13.706612275085302
Mean Absolute Percentage Error: 0.9731383925373656

test.loc[:, 'gdp_growth_pred'] = predictions

mape_by_country = {}
for country in countries:
    country_data = test[test[f'economy_{country}'] == True]
    mape_by_country[country] = mean_absolute_percentage_error(country_data['gdp_growth_lag'], country_data['gdp_growth_pred'])*100  # Convert to percentage

# Convert the dictionary to a DataFrame
mape_by_country_df = pd.DataFrame(list(mape_by_country.items()), columns=['Country', 'MAPE'])
mape_by_country_df = mape_by_country_df.sort_values(by='MAPE', ascending=True)

# Plotting MAPE by Country
plt.figure(figsize=(12, 8))
sns.barplot(x='MAPE', y='Country', data=mape_by_country_df, color='cornflowerblue', legend=False)
plt.title('Mean Absolute Percentage Error (MAPE) by Country')
plt.xlabel('MAPE')
plt.ylabel('Country')
plt.axvline(x=mape*100, color='red', linestyle='--', label='Overall MAPE')
plt.legend(loc='upper right')
for index, value in enumerate(mape_by_country_df['MAPE']):
    plt.text(value, index, f'{value:.2f}', va='center', ha='left', fontsize=9)

# Calculate RMSE by Country
rmse_by_country = {}
for country in countries:
    country_data = test[test[f'economy_{country}'] == True]
    rmse_by_country[country] = np.sqrt(mean_squared_error(country_data['gdp_growth_lag'], country_data['gdp_growth_pred']))

# Convert the dictionary to a DataFrame
rmse_by_country_df = pd.DataFrame(list(rmse_by_country.items()), columns=['Country', 'RMSE'])
rmse_by_country_df = rmse_by_country_df.sort_values(by='RMSE', ascending=True)

# Plotting RMSE by Country
plt.figure(figsize=(12, 8))
sns.barplot(x='RMSE', y='Country', data=rmse_by_country_df, color='springgreen')
plt.title('Root Mean Squared Error (RMSE) by Country')
plt.xlabel('RMSE')
plt.ylabel('Country')
plt.axvline(x=np.sqrt(mse), color='red', linestyle='--', label='Overall RMSE')
plt.legend(loc='upper right')
for index, value in enumerate(rmse_by_country_df['RMSE']):
    plt.text(value, index, f'{value:.2f}', va='center', ha='left', fontsize=9)

top_six_countries = mape_by_country_df.head(6)
bottom_six_countries = mape_by_country_df.tail(6)

# List of Top and Bottom Countries by MAPE
top_bottom_countries = pd.concat([top_six_countries, bottom_six_countries])['Country'].values.tolist()

# Generating Predcitions on the entire dataset for the top and bottom countries
predictions_all = best_model_proj.predict(data_lag.drop(columns=['gdp_growth_lag']))
full_data_preds = data_lag.copy()
full_data_preds['gdp_growth_pred'] = predictions_all

# Adding Country Names to the DataFrame
full_data_preds['Country'] = data_no_dummies['economy'].values

# Filter the DataFrame for the top and bottom countries
filtered_data = full_data_preds[full_data_preds['Country'].isin(top_bottom_countries)][['Country', 'gdp_growth_pred', 'gdp_growth_lag']]

# Creating Line Plot Grid for Top and Bottom Countries
plt.figure(figsize=(15, 15))
for i, country in enumerate(top_bottom_countries):
    plt.subplot(4, 3, i + 1)
    country_data = filtered_data[filtered_data['Country'] == country]
    plt.axhline(y=0, color='black', linestyle=(0, (5, 10)), linewidth=0.5)
    plt.plot(country_data.index, country_data['gdp_growth_pred'], label='Predicted GDP Growth', color='blue')
    plt.plot(country_data.index, country_data['gdp_growth_lag'], label='Actual GDP Growth', color='green')
    plt.fill_between(country_data.index, country_data['gdp_growth_pred'], country_data['gdp_growth_lag'], color='lightgray', alpha=0.5)
    plt.axvline(x=dt.datetime(2013, 1, 1), color='red', linestyle='--', label='Train/Test Split (2013)')
    plt.title(f'{country} GDP Growth Prediction ({"In Top 6" if country in top_six_countries["Country"].values.tolist() else "In Bottom 6"})')
    plt.xlabel('Year')
    plt.ylabel('GDP Growth Rate (%)')
    plt.legend()
    plt.xticks(rotation=45)
plt.tight_layout()
plt.suptitle('GDP Growth Rate Predictions for Top and Bottom 6 Countries by MAPE', fontsize=16, y=1.02)
plt.show()

Investing the Use of XGboost for GDP Growth Projection¶

Research Question:¶

Machine Learning Model Used: XGBoost¶

Our Dataset : World Bank¶

Dataset Overview¶

Reason for Selection¶

Features¶

Time Span:¶

Data source¶

Training The Model¶

Data Overview¶

Testing XGBoost Model With Basic Predictions¶

Baseline¶

Model Forcasting GDP Growth Rate Using Last Years Estimates¶

Hyperparameter Sensitivity Analysis (Max Depth)¶

Final Model¶

Breaking Down Models The Results¶

Mean Absolute Error (MAPE) & Root Mean Squared Error (RMSE) By Country¶

Comparing Prediction vs Actual GDP Growth Rate¶

Analysis of the Results¶

Challenges¶

Potential Future Questions:¶

Model Card¶

Model Card for GDP Growth Prediction using XGBoost (TEXT)¶

Model Overview¶

Model Description¶

Training Data¶

Evaluation Metrics¶

Limitations¶

Usage¶